home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / diffeq_45.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  220 lines
  1. ; $Id: diffeq_45.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ; Copyright (c) 1991-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. pro diffeq_45, funct, init, start, finish, times, yvalues, tol=tol,  $
  7.                report = report, Params = params, Listname = listname,$
  8.                Depvar = depvar
  9. ;+
  10. ; NAME:
  11. ;    DIFFEQ_45
  12. ;
  13. ; PURPOSE:
  14. ;    Solve system of first-order, ordinary differential equations:
  15. ;                    yi' = fi(t, y1(t),...yn(y)), i = 1,..., n
  16. ;                    ai  = yi(start),             i = 1,..., n         
  17. ;    using the Runge-Kutta method of order 4 and 5. Step size is selected
  18. ;    automatically and hence is variable.
  19. ;         
  20. ; CATEGORY:
  21. ;       Mathematical Functions, General
  22. ;
  23. ; CALLING SEQUENCE:
  24. ;       DIFFEQ_45, Funct, Init, Start, Finish, Times, Yvalues, $
  25. ;                  TOL = Tol, PARAMS = Params, REPORT = Report, $
  26. ;                  LISTNAME = Listname, DEPVAR = Depvar
  27. ;
  28. ; INPUTS:
  29. ;    Funct:    A character string containing the name of the user-supplied
  30. ;                  function implementing f = [f1, ...., fn].  This function
  31. ;                  should be written in IDL and have two arguments -- the scalar-
  32. ;                  valued time argument t, and the vector argument
  33. ;                  [y1(t), ... , yn(t)].  Additional constant parameters may be
  34. ;                  supplied through the keyword PARAMS.
  35. ;
  36. ;    Init:    The vector [a1, ..., an] of initial values.
  37. ;
  38. ;    Start:    The initial value of t.
  39. ;
  40. ;    Finish:    The final value of t. 
  41. ;  
  42. ;  
  43. ; KEYWORD PARAMETERS:
  44. ;    TOL:    The error tolerance.  The default is 1.e-6.
  45. ;
  46. ;    PARAMS:    A keyword to be passed to the function f. Params can be used
  47. ;        to specify constant-paramter values if f is a parametric
  48. ;        family of functions. See the example for the procedure 
  49. ;        DIFFEQ_23.
  50. ;
  51. ;        If the IDL function to compute f does not accept the keyword
  52. ;        PARAMS, then PARAMS should not be set in the call to DIFFEQ_45.
  53. ;
  54. ;    REPORT: If set, this flag signals that, at each step, the time
  55. ;        value, step size, and dependent variable values should
  56. ;        be written to the screen or to a file specified by keyword
  57. ;        LISTNAME.    
  58. ;
  59. ;     LISTNAME:    The name of the file to receive any output. The default is
  60. ;        to write to the screen. 
  61. ;
  62. ;    DEPVAR:    A string array of the names of the dependent variables to
  63. ;        be used in the output. Depvar(i) = name of variable i.
  64. ;
  65. ; OUTPUT PARAMETERS:
  66. ;          Times:    A vector of times at which f is computed.
  67. ;
  68. ;     Yvalues:    An array of y values. If ti = times(i),
  69. ;
  70. ;                 Yvalues(*, i) = f(ti,y1(ti),..., yn(ti))
  71. ;                               = [f1(ti,y1(ti),..., yn(ti)),  ..., 
  72. ;                                              fn(ti,y1(ti), ...,yn(ti))].
  73. ;
  74. ; COMMON BLOCKS:
  75. ;    None.
  76. ;
  77. ; SIDE EFFECTS:
  78. ;    None.
  79. ;
  80. ; RESTRICTIONS:
  81. ;    None.
  82. ;
  83. ; EXAMPLE: 
  84. ;       See the example for DIFFEQ_23.
  85. ;
  86. ; PROCEDURE:
  87. ;           The constants of Fehlberg are used to obtain the Runge-Kutta
  88. ;           formula. See "Klassiche Runge-Kutta Formeln vierter und
  89. ;           niedrigerer Ordnung mit Schrittweiten-Kontrolle und ihre
  90. ;           Anwendung auf Warmeleitungsprobleme", Computing 6, 1970,
  91. ;           pp. 61 - 71.
  92. ;
  93. ; MODIFICATION HISTORY:
  94. ;                CAB, Sept., 1991.
  95. ;-              
  96.  
  97. On_Error,2
  98.  
  99. ; Check parameters
  100. if n_params(0) lt 4 THEN $
  101.    message, "Missing parameters"
  102.  
  103.  
  104. if KEYWORD_SET(listname) THEN $
  105.    openw, unit, /get, listname $
  106. else unit = -1
  107.  
  108.  
  109. if KEYWORD_SET(tol) eq 0 THEN tol = 1.e-6
  110. if KEYWORD_SET(Params) eq 0 THEN  $
  111.   Paramset = 0   $
  112. else Paramset = 1
  113.  
  114.  
  115.  
  116. if KEYWORD_SET(report) eq 0 THEN report = 0 $
  117. ELSE  BEGIN
  118.    SN = size( depvar)
  119.    n = n_elements(init)
  120.    if (SN(1) eq 0) THEN BEGIN
  121.       I = indgen(n)
  122.       Names = ['Var' + StrTrim(I, 2)]
  123.    ENDIF ELSE  $
  124.       if SN(1) lt n THEN BEGIN
  125.       I = Indgen(N)
  126.       Names = [DepVar, 'Var' + StrTrim(I(SN(1) : N-1),2)]
  127.       ENDIF else Names = depvar
  128.   
  129.    printf, unit, format = '(A13, 2x, A13, 2x, $)', "Times", "Stepsize"
  130.    for i = 0, n-2 do printf,unit,format = '(A13,2x,$)',Names(i)
  131.    printf,unit, format = '(A13)',Names(n-1)
  132.   Printf,unit, " "
  133. ENDELSE
  134.  
  135. ; Constants of Fehlberg
  136. a = [1/4., 3/8., 12/13., 1., 1/2.]
  137.  
  138. b = [ [ 1./4, 0, 0, 0, 0],  $
  139.       [ 3./32., 9/32., 0, 0, 0], $
  140.       [1932./2197, -7200/2197., 7296/2197., 0, 0], $
  141.       [439/216., -8., 3680/513., -845./4104., 0], $
  142.       [-8/27., 2., -3544/2565., 1859./4104, $
  143.                     -11/40.]]
  144.  
  145.  slope_weight =   [16./135., 0, 6656./12825., 28561/56430., -9./50., 2./55.]
  146.  err_weight = [-1./360., 0, 128./4275., 2197./75240., -1./50., $
  147.                  -2./55.]
  148.  
  149. ; Initialize
  150. h = (finish - start)
  151. minh  = h/20000.
  152. maxh  = h/5.
  153. times = start
  154.  
  155. h = h/100.
  156. t = start
  157. v = init
  158. yvalues = v
  159.  
  160. if report  ne 0 THEN BEGIN
  161.    printf,unit, format ='(G13.6, 2x, G13.6, 2x, $)', t, h
  162.    for i = 0, n-2 do printf,unit, format ='(G13.6, 2x,$)', v(i)
  163.    printf,unit,format = '(G13.6)', v(n-1)
  164. ENDIF
  165.  
  166.  
  167. slopes = fltarr(N_ELEMENTS(v), 6)
  168.  
  169. errbound = tol * max([sqrt(total(v^2)), 1])
  170.  
  171. ; compute yvalues for smaller (as needed) step sizes
  172.  
  173. while t lt finish and h ge minh DO BEGIN
  174.  
  175.  if t+h gt finish THEN h = finish - t
  176.  
  177.  if Paramset eq 0 THEN BEGIN
  178.    slopes(0, 0) = CALL_FUNCTION( funct, t, v)
  179.    for i = 1,5 do $
  180.        slopes(0, i) = CALL_FUNCTION( funct, t + a(i-1) * h, v +   $
  181.                                     h * (slopes(*,0:4) # b(*,i-1)))
  182.  ENDIF ELSE BEGIN
  183.    slopes(0, 0) = CALL_FUNCTION( funct, t, v, Params = Params)
  184.    for i = 1,5 do $
  185.        slopes(0, i) = CALL_FUNCTION( funct, t + a(i-1) * h, v +   $
  186.                                     h * (slopes(*, 0:4 ) # $
  187.                                     b(*,i-1)), Params =Params)
  188.  ENDELSE
  189.  
  190.  err = ( h * (slopes # err_weight))^2
  191.  
  192.  err = sqrt(total(err))
  193.  err_bound = tol * max([sqrt(total(v^2)), 1])
  194.  
  195.  if err le err_bound THEN BEGIN
  196.     t = t + h
  197.     v = v + h*(slopes # slope_weight)
  198.     times   = [times,t]
  199.     yvalues = [[yvalues], [v]]
  200.     if report  ne 0 THEN BEGIN
  201.       printf,unit, format ='(G13.6, 2x, G13.6, 2x, $)', t, h
  202.       for i = 0, n-2 do printf,unit, format ='(G13.6, 2x,$)', v(i)
  203.       printf,unit,format = '(G13.6)', v(n-1)
  204.    ENDIF
  205.  ENDIF
  206.  
  207.  if KEYWORD_SET(do_print) ne 0 THEN $
  208.     print, t, h, v
  209.  
  210.  if err ne 0 THEN  $
  211.     h = min([maxh, .8*h*(err_bound / err)^(1/5.0)])
  212.  
  213.  endwhile
  214.  
  215.  if ( t lt finish) THEN $
  216.      printf,unit, " Beware of singularity"
  217.  if unit ne -1 THEN Free_lun,unit
  218. return
  219. end
  220.